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This  paper  describes  algorithms  to  elieit  multinomial  sampling  on  a 
computer  in  ways  that  protect  the  sampler  against  an  excess  computing  cost 
per  sample.  Section  2 presents  the  multinomial  model  together  with  an  equiv- 
alent representation  in  terms  of  a series  of  binomial  sampling  experiments. 

The  binomial  representation  is  further  discussed  in  Section  4.  Sections  3 
and  5 demonstrate  the  danger  of  using  too  simplistic  a sampling  scheme  if 
execution  time  is  a concern.  Section  6 describes  how  a normal  approximation 
to  the  binomial  distribution  can  make  execution  time  virtually  independent 
of  n.  A criterion  of  acceptability  is  described. 

Section  7 describes  .in  acceptance-reject  ion  technique  that,  when  used 
with  a Poisson  sample,  allows  the  desired  binomial  sampling  exactly.  Using 
.1  normal  approximation  to  the  Poisson  distribution,  oni-  can  again  gener.ite 
multinomial  samples  virtually  independent  ol  n.  A criterion  of  acceptability 
is  described.  The  method  of  Section  7 may  apply  with  given  accuracy  in  c.isi's 
in  which  the  normal  approximation  to  the  binomial  distribution  does  not  apply. 

Section  8 describes  a procedure  for  binomial  sampling  based  on  the  inverse 
transform  method.  Although  the  mean  execution  time  is  proportional  to  n,  the 
procedure  is  intended  for  small  n.  Section  9 describes  the  ordering  of  the 
serial  binomial  experiments  in  response  to  alternative  objectives.  Section  10 
describes  algorithm  Ml  which  puts  all  the  suggestions  of  earlier  sections 
together.  The  section  also  describes  procedures  for  reducing  the  expense  of 
testing  to  see  if  cost-saving  generation  methods  apply. 
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Introduction 

Multinomial  sampling  Is  an  Inherent  feature  of  many  computer 
based  simulations.  The  most  notable  types  of  simulation  that  rely 
on  this  form  of  sampling  have  to  do  with  population  growth.  For 
exain)lc, given  n Individuals  In  an  Initial  state,  one  may  wish  to 
determine  the  numbers  that  pass  to  new  states  l,...,k  in  a unit  time 
period.  Here  p^^,  the  probability  of  passing  to  state  i, 
may  be  constant  for  the  tenure  of  the  simulation  run  or  may 
change  in  each  unit  time  period  In  response  to  environmental  change. 

Let  denote  the  number  that  pass  to  state  1.  Then  in  a given 

year  with  specified  n and  p^,...,Pj^,  Xj,...,Xj^  form  a multinomial  sanipU'. 

Although  this  form  of  sampling  is  hardly  rare,  little  If  anv  discussion 
has  appeared  In  the  published  literature  regarding  i'ow  to  perform  tiiis 
sampling  efficiently  on  a computer.  This  paper  addresses  the  problem 
of  computational  efficiency  and  describes  algorltliras  to  effect  multinomial 
sampling  In  ways  that  protect  the  sampler  against  an  excess  cost  per 
sample . 

Section  2 presents  the  multinomial  model  together  with  an  equivalent 
representation  In  terms  of  a series  of  binomial  samplittg  I'xper Inu'iit s . 

The  binomial  representation  is  further  discussed  in  Section  A.  Sections 
3 and  5 demonstrate  the  danger  of  using  too  simplistic  a sampling  scheme 
if  execution  time  Is  a concern.  Section  b describes  how  a normal  approx- 
imation to  the  binomial  distribution  can  make  execution  time  virtually 
independent  of  n.  A criterion  is  presented  there  for  determining  when 


the  approximation  applies  with  a specified  upper  bound  on  absolute  error. 
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The  criterion  Is  In  terms  of  an  inequality  that  easily  can  be  built 
into  a computer  program. 

Section  7 describes  an  acceptance-rejection  technique  that,  when 
t used  with  a Poisson  sample,  allows  the  desired  binomial  sampling  exactly. 

I Using  a normal  approximation  to  the  Poisson  distribution,  one  can  again 

generate  multinomial  samples  virtually  Independent  of  n.  A criterion 
of  acceptability  is  described.  The  method  of  Section  7 may  apply  with 
given  accuracy  in  cases  in  which  the  direct  approximation  of  the  binomial 
distribution  does  not  apply. 

Section  8 describes  a procedure  for  binomial  sampling  based  on  the 
inverse  transform  method.  Although  the  mean  execution  time  is  pro- 
portional to  n,  the  procedure  is  Intended  for  small  n for  which  evidence 
in  [ 4 1 indicates  its  appeal.  Section  9 describes  the  ordering  of 

! the  serial  binomial  experiments  in  response  to  alternative  objectives. 

For  example,  if  one  samples  in  order  of  decreasing  p^,  the  odds  that 
the  normal  approximation  in  Section  6 applies  are  enhanced.  Section  10 
* describes  algorithm  M3  which  puts  all  the  suggestions  of  earlier  sections 

together.  The  section  also  describes  procedures  for  reducing  the  expense 
I of  testing  to  see  if  cost-saving  generation  methods  apply. 


2.  The  Multinomial  Model 

Consider  a series  of  n Independent  trials  on  each  of  which  just 
one  of  k mutually  independent  events  , . . . , Ej^  occurs.  The  event 
F.j  occurs  with  probability  p^.  Let  denote  the  number  of  time  that 

occurs.  Then  the  joint  probability  mass  function  (p.m.f.)  of 
X - (Xj X^)  Is 

A 
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( 1 ) 


pr  (X^  » ; 1 


k 

1, . . . , k)  - n!  77  (P,/"! ’) 

i-1  ^ 


0 < p^,  0 5 Xj  i 


1 


n. 


Here  X has  the  multinomial  distribution  denoted  by  M(n,  p^ , . . . , p^^). 
As  an  example,  consider  a population  of  Immature  female  elephants  of 
the  same  age  in  a park  preserve  in  Africa.  During  a given  year  each 
elephant  experiences  one  of  k ■ 4 events; 


“ dies. 

= survives  and  remains  immature. 

E^  “ survives  and  matures. 

E^  - survives,  matures  and  conceives. 

Maturity  means  capable  of  conceiving.  Since  these  categories  are 
exhaustive  p^^  + P2  + p^  + p^  * 1. 

Let  Ej^  denote  the  occurrence  of  an  event  other  than  E^  and  let 
k 

M - I X 1 - 1,...,  k 

j-1  J 
Jl‘1 

denote  the  number  of  trials  on  which  E^  occurs.  Here  E^  has 
probability  ^tid  X^  and  have  the  joint  p.m.f. 

( 2)  pr(X^  - X,  - n-x)  - (")  p^d-p^)"'*  x - 0,...,  n. 

But  (2)  Indicates  that  X^  has  the  binomial  distribution  B(n.  P|)- 

Suppose  that  events  , . . . , E^  occur  Xj^,...,Xj  times,  respectively. 


for 


j < k and  define 
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I 

I 

1 

I 
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(3)  Z,  - I ^ X • 

J i-1  ^ 

Then  one  can  show  that 


shortly  these  results  play  central  roles  In  multinomial  sampling. 

3.  Simple  Random  Sampling 

Let  us  now  recast  the  problem  In  the  more  familiar  terms  encountoroci 
in  discrete  event  simulation.  Suppose  that  a given  state  has  a pop- 
ulation size  n and  that  each  member  has  a probability  p^  of  going 
to  state  1 for  1 “ 1,...,  k.  The  problem  then  becomes  one  of 
determining  X^,  the  number  of  Individuals  that  move  from  the  given 
state  to  state  1.  More  generally  one  wants  to  determine 
X " (Xj  > • • . • Xj^)  . 

Many  methods  come  to  mind  for  sampling  X from  M(n,  p^,...,  p^) 

on  a computer.  The  simplest  relies  on  Independent  random  sampling  of 

the  n events.  Let  be  a uniform  deviate  on  (0,1)  whose 

distribution  Is  denoted  by  U(0,1).  Then  algorithm  Ml  affects  the 

desired  sampling.  Here  I,  , (U.)  denotes 

(qj.l.qjl  1 


1 

I 

j 
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i 
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till"  iiuiicator  tuiu’tlou 


^ ^ \ n 1 

(q  . q ) 1 


r 


<I  ' U • q 

*1-1  i i 


lit  liorw  i si’  . 


Algol  it  lim  Ml 


1. 

Xj  ' 0 

1 = 1 k. 

T 

1 - 1. 

3. 

Sample  l!^ 

from  11(0,1). 

4. 

” • '’'(.li 

1 

i 

5. 

X - X + 1 
m m 

• 

6. 

If  i “ n 

deliver  X = (X^ 

* • • • 

, x^). 

7. 

i « i + 1 . 

8. 

Go  to  3. 

l.et  T(B) 

denote  the 

mean  execution 

t ime 

of  algorithm 

mean  time 

to  sample  X 

from  Ml  has 

the 

form 

T(>!I)  * + nf(k,p)a^  wtioro  aj  and  a^  depond  on  computer  and 

programming  considerations  and  t'(k,p)  depends  on  k and 
p = (Pp...,  Pj^)  tlirougii  the  algoritlim  selected  lor  executing  step  4, 
I’rocedures  TJ  and  T4  in  Fishman  1 5,  Sec.  9.18]  describe  tabjjng 

algorltlims  based  on  MarsagHa  [ 9)  that  make  Ilk,  p)  independent  ol  k. 
These  procedures  are  tiighly  efficient  timewise  but  require  more  space 
tlvin  search  algorithms  such  as  procedures  T1  and  T2  in  | 5.  Fee.  d.lg] 

do.  Since  all  these  algorithms  for  executing  step  4 iiave  their  greatest 
appeal  when  is  fixed  throughout  a simulation,  a question  remains 

reg.irdlng  their  suitability  for  a simulation  in  which  state  transition 
probabilities  change  as  time  elapses.  More  importantly  the  linear 
growth  of  T(Mn  with  n makes  Ml  prohibi t ivel v expensive  for  large  n. 
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4.  Sequential  Binomial  Sampling 

Among  the  alternatives  to  simple  random  sampling,  most  rely  on 
some  form  of  binomial  sampling.  Algorithm  M2  offers  a prototype. 

Let  T(M2,  B)  denote  the  mean  execution  time  for  M2  using  j 


Algorithm  M2  (n,£) 


1. 

Z 

0. 

2. 

"l 

*■  0 for 

1 * ka 

3. 

1 

1. 

4. 

Sample  from  8(n-Z,  Pj/(l-q^  j^)) 

5. 

If 

= n-Z, 

deliver  X^^ , . . . , X^^. 

6. 

Z « 

Z + X^. 

7. 

If 

i = k-1. 

X^  n - Z and  deliver 

8. 

i t 

1+1. 

9. 

Go 

to  4. 

Algorithm  B (to  be  described)  for  sampling  from  f^CN,  p) . Then 


k-1 

( 6 ) T(M2,  B)  = a + a.k  + a,  I Ell)(n-Z  , pj.  B) ) . 

1 J A 3 j^l 

I 

1 ^0“  ° Pi  * ^ ^ 

vdiere  E[d(N,  p,  B)J  is  the  mean  sampling  time  for  algorithm  B given 
N and  p,  Is  defined  In  (3  ),  E Is  the  expectation  operator 


and  a^,  a^  and  a^  depend  on  computer  and  programming  considerations. 
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5.  Bernoulli  Tr lals 

Notice  that  ( f< ) always  has  a term  proportional  to  k. 

However,  the  form  of  dependence  on  n depend  on  the  selected  binomial 
sampling  algorithm.  Ahrens  and  flieter  [1,2]  and  Fishman 
[ 3 , Sec.  9.14]  describe  several  alternative  algorithms.  If  one  uses 
independent  Bernoulli  trials  to  sample  X from  B(N,  p)  as  in 
algorithm  BE,  one  has 


r>(N,  p,  BE)  ' min  (p,  1-p)  + b^N, 


Algorithm  BE  (N , p) 


1. 

X 0. 

2. 

P ^ P • 

3. 

If  p > 

0.5 

p *■  1-p. 

4. 

i - 1. 

3. 

Sample 

U from  t/(0,l). 

6. 

If  U < 

P. 

X X+1. 

7. 

If  1 < 

n. 

1 *■  1+1  and 

8. 

If  p < 

0.5 

deliver  X. 

9. 

X - n-X. 

10. 

Pe liver 

X. 

t ies 

*’2 

and 

bj  depend  o 

considerations.  The  quantity  b^  denotes  the  time  spent  executing  a 
trial,  regardless  of  its  outcome,  whereas  b^  denotes  the  time  spent 
incrementing  X when  a success  occurs.  Generally  b^  >>  b^.  Steps 
2,  3,  8 and  9 also  require  explanation.  One  can  easily  show  that  If 
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r ^ 

1 

! 

i 

j 

1 

X Is  from  B(N,  p)  then  N-X  is  from  B(N,  1-p) . By  usfng 
p " min(p,  1-p),  the  mean  number  of  additions  In  step  5 Is 
Np  •s  Np . For  example,  p = 0.99  gives  Np  = N ''  0.1  versus 
Np  " N X 0.99.  Although  the  cost  of  addition  Is  snuall,  nevertheless 
the  saving  Is  worth  making. 

If  one  follows  M2  then  the  mean  time  to  sample  X^  is 
(8)  F.lD(n-Z^_j,  p^,  BE)]  = ^2"  '“9^)  + 

The  possibility  of  not  executing  step  4 because  j ^ contributes 


the  coefficient  1 

P 

-‘>1- 

^ to  bj^.  One  now  can  show 

that 

k-2 

n 

‘'i 

( 9 ) T(M2,  BE)  = 

“3  - 

^ ‘’l  " ^ ^ 

V 

L 

1 = 1 

k-1 

k-1 

+ n a^  (b^ 

1 

i=l 

mln(p^,  l-q^)  + b^(k-l)pj^  + 

^ ) 

' 1 p,  ] • 

'1-1 

Since  the  term  in 

n 

dominates,  the  Bernoulli  method  also  loses  appeal 

as  n increases. 

6.  Normal  Approximation 

Provided  that 

n 

is  sufficiently  large  for  a 

given 

p,  one  can 

turn  the  large  n to  a virtue  in  multinomial  sampling.  It  is  a well 
known  result  in  statistics  that  if  X is  from  BCN,  p)  then 
Y = (X-Np) /*4lp(l-p)  has  a distribution  that  converges  to  the  normal 

distribution  V(0,1)  as  N increa.ses.  U-t  Y be  fromN(0,l).  T''-'n 
one  computes  X as 

(10)  X • max  [0,  min(N,  <Y/Np(l-p)  + Np  + 0.5>)]  . 

Algorithm  BN  describes  how  to  effect  this  sampling. 

^The  quantity  <0>  denotes  the  largest  integer  in  0. 
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Al>;ori  t lim  BN  ( N , p ) 

1 . a - Np . 

.>.  h < ..i('l-p')  . 

i.  a a + 0.05. 

4.  Sample  V from  V(0,1). 

5.  X max[0,  min(i'J,<  Vb  + a >)]. 

6.  Deliver  X. 


For  p the  t-rror  in  treating  as  binoniiai  dei' riMSi-s  as 

N increases.  Assi’ssnient  ol  tin's  erroi  r.ees  as  lar  hack  as  lls|>ensk\’ 

[10] .  Here  we  use  a result  in  Makabe  [/]  that  provides  a tighter  erroi 
bound  than  Uspensky  does.  Let 

(11)  X - 1 = y/^pq  + Np-0.5  q=l-p,  P<q. 


Then  for  Npq  2 25 

(12)  i 1 pr  (X  S x-1)  - pr  (Y  < y) | 


where 

q-p  , 0.053  + 0.027(p-q)  + 0.054(p-q)^  , ^-1.5 

^ ' A 'nr- 

o>'2TTNpq 

Table  1 shows  the  minimal  N required  to  satisfy  (13) 
for  .Ai  = 0.01  and  0.02.  The  results  were  obtained  by  appl  icat  ion  oi 
the  Newton-Raphson  method.  Wi'  discuss  the  usi>  ol  (1  i)  in  prai  l lie  in  So<  I ion 


10, 

For  an  X generated  by  alogrithm  BN’,  one  can  slu  w that 
the  mean  execut  ion  time-  of  BN  is 
(14)  D(N,p,BN)  = b^. 

If  BN  is  used  in  step  4 of  M2,  one  has 

T(M2.  BN)  = a^  + a^k  + a,^b^(k-l). 


(15) 


since  a little  work  shows  that 


(Ih)  llm  T(M2,BN)  = ^>3  + a^k  + a^b^^), 

n »«' 

clearly  one  can  effevt  multinomial  sami>linn  lelativi'lv  i miepeiuient  ol 
n,  proviileii  that  n is  siif  f Ic  lent  1 v larpe.  Ttie  kev  componiMit  s ol 
tills  result  are  twofold.  Klrstly,  tlie  converpa’iice  of  tlie  hinomial  to 
tile  noriOi'il  distribution  is  a necessary  inp.ri'd  ient  . Secondly,  t ■e 
independence  of  the  cost  of  normal  sampllnp  of  the  binomial  parami'ters 
in  step  4 of  BN  is  critical. 

If  one  uses  BN  exclvisivelv  in  M2,  a serious  error  can  arisi' 

In  practice.  Since  successive  calls  to  BN  use  ii,  n-i:|,  n-z,,...,  n-;’.|^_j 

it  is  entirely  conceivable  that  (n,  pi),  (n-z, , p'  (u-o  p’  ) lot 

i 2 ‘ i ’ i + 1 

i < k-1  satisfy  (12)  for  a given  A but  (n-z,.,,  p 1+ >),...,  (n-.:  p'  ) 

1+1  , I j 

ilo  not.  When  tlie  constraint  on  n-z  Is  lu’t  mi-t,  oiu'  lu'ods  .111  .1lien1.11  ivi' 

i 

sampling  procedure  to  di'termine  X , . . . Tlieretore,  oiu'  laniuU  re  1 v I'li 

BN  to  solve  all  needs.  We  describe  a 1 1 1'rn.it  i ves  In  Sect  iiMis  7 and  S. 


7 . 1’ o I s son  No Jt  ho d 

'I'ac  problem  one  now  faces  is  to  find  an  elficient  wav  ol  s.impling 
from  B(N,  p)  when  N does  not  satlsfv  a specilied  A|  in  fl.’)  li'i  .1 
given  p.  Ahreii.s  and  Dieter  ill  .ind  12|  .uni  Klshm.in  |S,  See.  S.l.),  and 
[ 4 1 describe  alternative  methods  of  binomial  sampling.  Among  these 
tile  I’ols^on  method  in  I 4 1 and  [ ' | seems  most  appealing  tor  large 
N wtien  BN  d.  es  not  applv.  let  X be  a I’olsson  distributed  random 


variable  with  p.m.f. 


( I /)  pr(X  - x) 


- 11  X , , 

e |i  /xl 


X ^ 1) , I . . . . 


I't  lii'i  w 1 se  , 


denoted  bv  PCp).  l.et  II  be  from  (1(0,  1).  It 


X <;  N 


.N-X-<X-* 


(N-x): 


\ li(l-p)/p 


then  X Is  from  i^N,  p) . Since 


pr(S2|Sj)  “ /N'.e’'  (l-p)^X  ' pr(S|>, 


the  mean  number  of  X's  from  (17)  need  to  find  an  X from  8(N,  p) 


(18)  ^ ° l/pr(S,,  |Sj)pr(Sj)  ” Nl  e‘'(l-p)’^X''“N/vA>! 


This  expression  Is  minimized  by 


N - < N(l-p)  > 


N(l-p)--.N(l-p)  V p 


P(  «-N(l-p)  - +l)/(l-p) 


otherwise. 


See  t"’  • Sec.  9.141  or  [ 4 ]. 

Now  the  properties  of  (18)  are  of  seminal  importance.  First 1v, 


1 Im  C (\i,  p)  - 1 / (1-p) 


u 


Secondly,  Table  2 shows  p)  (I  - fer  selected  N 

and  p.  The  results  reveals  a virtual  insensitivity  of 

1 /2 

C^(p*,  p)  to  N and  that  p)  ~ l/(l-p)  • Thirdly,  if 

one  makes  use  of  the  observation  that  N - X is  from  B(N,  1-p) 
wlten  X is  from  B(N,  p),  then  one  can  replace  p bv 
p » min(p,  1-p)  everywhere  in  this  section  and  use  X if  p • O.S 
and  N - X if  p > 0.5.  This  third  property  allows 

(21)  ~ l/(l-p)^^^  *'  ~ 1.41, 

where  p replaces  p in  (18). 

Now  the  remaining  cost  of  the  Poisson  method  resides  in  the 

algorithm  selected  for  Poisson  sampling.  Alirens  and  Dieter  [l  1, 

and  [ 2 ] and  Fishman  ( A ] and  ( 5 , Sec.  9.1JI  describe  algoritlims 

1/2 

whose  mean  execution  times  are  proportional  to  (u*)  . Since 

VI*  = N - <N(l-p)^  ~ Np  for  large  N,  one  can  at  least  reduce 

the  cost  of  binomial  sampling  and  consequently  of  multinomial 

1/2 

sampling  to  a factor  proportional  to  N . This  clearly  improves 
on  the  linear  dependence  on  N in  random  and  Bernoulli  sampling. 

In  practice,  one  also  can  Improve  on  these  square  roiH  methods 
considerably  by  exploiting  the  limiting  behavior  I'f  P (vi ) . It  X is 
from  P(p),  then  (X-vi)//vr  has  a distribution  that  converges  to 
N/(0,1)  as  VI  Increases,  The  error  or  approximation  has  been  stndieil 
by  a number  of  writers,  Including  Cheng  [3J  and  Makabe  and  Morlmura 
[s].  Here  we  use  the  result  of  Miikabe  and  Morlmura,  who  give  the  tigliter 
bound.  Let  x be  an  integer  and  let  y satisfy 

(22)  X - y + ii  -r  0.5. 


Let  Y be  from  W(0,l)i  Then  for  V»  1 


IS 


(21) 


WluTf 


^2  ''  |pr(X  • x-1 ) - pr  (Y  < y)  | 


- + O.OSAA/ii  ♦ l).lMOS/ii  ().2/A)/ii' 

6 ►' JjTIp 

+ 0.0065/v.^'^^  + (1  + *.“2  Ai 


Wf  return  to  thio  issue  et  spi-eitvin>;  .in  i.i'pei  lunind  on  I’lior 

!ii  Seet  ion  10^  when  bringing  the  alternative  binomial  algorithms  iiigi-ther 
in  an  omriibns  multinomial  s.impling  algorithm. 

Algorithm  Bl’N  d.*.si-r ibt-s  luiw  to  effeet  binomial  sampling  via 
the  I’oisson  met  liod  using  the  normal  approxinvit  ion . Here  step  19  tests 
S,  and  needs  explanation.  By  working  with  V » - In  U,  whieh 


is  an  exponential  deviate,  one  transforms 


S 


2 


to 


(23)  V - -In  U i ('\>  - N + X)  In  1 + In  1(N-X).'J  - In  (<)>.’)  , 


trom  which  step  19  follows.  Also,  by  using  a table  of  1+1  entrei's 
for  ln(i.')  tor  i L one  improvi-s  .■  I t u i.ncy . Ilu-  el.oi.e  of  i . ,s  i h.-  u.s,  i 
Steps  10  and  16  account  for  situations  in  whicli  the  tables  .ire  inade(|u.it  i'  by 
using  Stirllng'e  approx  im.it  ion  for  large  factorials.  Apart  from  the  time 
spent  executing  steps  10  and  lb  wlien  their  conditions  are  true,  .ilgorithm 
Bl’N  has  mean  execut  ion  t imv 


i 


1 


ib 


Algorithm  BPN  (N.  p) 

1 

Table:  f i 0,  f - In  j for  I - 1 1. 


c ' / 2n  . 

1. 

P 

P- 

2. 

If 

p > 0.5,  p ♦ 1-p  . 

3. 

r ♦ 

N(1  - p). 

4, 

s 

<r> . 

5. 

p * ♦ 

N-s. 

6. 

If 

r-8  > p,  p*  ♦ p(8+l)/(l-p). 

7. 

X -t- 

(1/p  - Dm*. 

8. 

♦ ^ 

In  X . 

9. 

m 

A 

V 

10. 

If 

m > I-l,  q ♦ c - m + (ro  + 0. 5)  In  m and  go  to  12. 

11. 

R 

^nrt-l  * 

12. 

a ■«- 

CM 

13. 

Sample  Y from  N(0,1). 

14. 

X 

max(0,<Ya  + p*  + 0.5>). 

15. 

W 

N - X. 

16. 

If 

W>I-1,  h'«-c-W+(W  + 0.5)*  InW  and  go  to  18 

17. 

h ♦- 

W. 

18. 

Sample  V from  E(l). 

19. 

If 

V > (m  - W)  (^  - g + h,  go  to  13. 

20. 

If 

p i 0.5  deliver  X. 

i 

1 


21.  Deliver  W. 


17 


(24)  K[D(N,  p,  BPN)1  - ~ 5 + /n>ax(p,  1-p) 


Tlien  M2  with  n,  Pj^,...,  p^^  and  using  BPN  has  mean  time 


r k-i 

(25)  T(M2,  BPN)  ~ + a^k  + A b^  (1  - q"_^) 


k-1 


+ ,)/n'ax(p  , 1 - q )1 

i = l 


1/2 


1 


^ a.  + a k + a (b  + /l  b )(k-l)  . 
3 4 5 5 6 


8.  Inverse  Transform  Method 

I 

There  remains  the  issue  of  how  to  sample  from  the  binomial 
distribution  when  neither  the  normal  approximation  of  Section  6 nor 
the  Poisson  method  of  Section  7 applies.  Timing  experiments  in  I 4 I 
revealed  that  the  inverse  transform  sampling  method  deserved  serious 
consideration  when  Np  < 15.  In  fact,  it  was  the  least  time  consuming 
for  most  (N,  p)  combinations  examined  wlien  compared  with  three 


I alternative  algorithms.  Although  one  may  justlfiablv  consider  using 

i 


L«*t  U be  a random  variable  from  U(0,1).  Since 

t 


. 1'  > 


X * 0,1..., 


one  can  determine  X a.s 


\“  min  (x  : 0 ^ t^j  x = 0,1,...)  . 


for  X from  8(N,  p)  one  has 


/ N - X 

**0  ■ ■ X + 1 


_i!L_ 

1 - p 


Algorithm  B1  [4]  describes  the  steps  need  to  sample  X from 

8(N,  p)  using  the  Inverse  transform  method.  Here  mean  execution  time 


Algorithm  B1  (N, 


A, 

B,  C, 

D are  double  precis 

1. 

q 

P- 

2. 

If 

p>0.5  q'^l-p. 

3. 

s 

1 - q. 

4. 

A ^ 

1. 

3. 

B .e 

q/s. 

6. 

C -e 

(N  + 1)B. 

7. 

D 

A. 

8. 

X 

0. 

9. 

Sample  U from  U(0,  1) 

10. 

V -e 

U/a  . 

11. 

If 

V 5 A go  to  16. 

12.  X ♦ X + I. 


13.  D t D(C/X  - B), 


lA. 

A - A + 

D. 

15. 

If  X < 

N 

16. 

If  p > 

0.5 

17. 

Del iver 

X. 

18. 

Deliver 

N • 

has 

the  form 

i2h) 

ti(N.  p, 

Bl>  = b^ 

+ b 

Then 

M2  wit  h 

n , P • * 

t(m: 

2,  Bl)  = a3 

+ a , k + 

*4 

«5 

;)iu1  using  bL  l';is  moan  exooution  limo 


[ i=l  ^ ‘ 1-1  ’ ' j 


•■I,  Ordering  the  StaJ:es  for  _Kxei'iit_ion 

As  the  virevious  sections  show,  viewing  multinomial  sampling  as  .i 
sequence  of  Binomial  sampling  experiments  enables  one  to  take  advaniai’i' 
of  BN,  wlien  it  applies,  then  Bl’N  wlien  it  applies  and  finally  Bl 
when  N mln(p,  1 - p)  is  relatively  small.  In  practice,  one  wants  to  induce 
additional  efficiency  by  thought fullv  choosing  the  order  in  which  states 
1 « 1 k are  sampled.  Let 


3 “ (J  j , • • • * .1  j^) 


0 ' .1^  •'  k 


JO 


denote  a vector  of  Integers  such  that  i j **  J £ ^ ^ 

1,  « “ 1,...,  W.  Here  each  can  assume  values  I through  k hut 

no  two  Jj'**  I’Moal.  The  objeitive  now  is  to  pick  in  an  optimal 

way  with  regard  to  a specified  norm. 

Let  us  first  concentrate  on  (13).  Clearly,  tite  larger  Npq  is, 
the  more  likely  one  is  to  satisfy  a specified  A.  Prior  to  generating 
the  sample  X ,X  ,...,X  , X has  variance  Np,  (1-p,  ) . One  irame- 


Jl  Jj Jk’  Ji 


.1 


Ji' 


diately  desirable  objective  is  to  choose  £ so  that 


This  follows  from  the  assignment 


i “ 1 k - 1. 


(J4) 


'’.i 


i = 1 k - 1, 


’ i ' 1+1 

The  choice  in  (29)  Is  also  compatible  with  m.iximlzing  tlie  odds  ot 
satisfying  (22)  for  a specified  A. 

If  neither  BN  or  BPN  apply,  then  one  needs  to  reconsider  the 
arrangement  J.  Clearly  the  rule 

p,  5 p 1 = I, . . . , k 


(30) 


P.  i P, 
Jk  '1 


i - 1 k - 2 


'l+l 


improves  on  (29)  by  assigning  the  potentially  most  costly  sampling  lo  i hi' 


kth  sampling  position.  Using  the  identity  a.  - rt  - X,  - X X 

■'k  'l  Jj  J 


k-1 


avoids  this  cost. 


10.  Putting;  It  All  To»;fttu‘r 


Since  alternative  forms  of  binomial  sampling  appear  attractive  for 
different  (N , p)  pairs,  one  would  like  a multinomial  sampling  procedure 
for  each  (N , p)  combination  encountered.  Algorithm  M3  does  this. 

Here  A*  and  A"  are  error  upper  bounds  on  (13)  and  (3d)  respec- 
tively. The  reader  should  note  that  BN,  BNP  and  B1  are  called  separati’lv 
rather  than  as  components  of  one  omnibus  binomial  procedure.  This  lorm 
is  intentional.  It  reduces  the  amount  of  testing  for  choosing  among  pro- 
i«'(.lu  r«’S . 


Algorithm  M3  (N,  p.  A',  A") 


(■iven  c = 1 /6  > 2 n . 


1.  Determine  1^,...,  such 


1 “ 

1.. 

..,  k - 1, 

2 . 

"i  ' 

0 

for  i = 

1 k. 

3. 

i ' 

1. 

4. 

q ' 

0. 

5. 

p 

"Ji 

6. 

N - 

n . 

7. 

A ' 

8. 

8. 

m * 

»Np(l  - p) . 

9. 

If 

m < 

5 go  to 

12. 

10. 

If 

sat isf ies 

cond it  ion 

BN(N 

. p) 

and  go 

to  21. 

11. 

A - 

11. 

12. 

m ♦ 

*Np . 

that 


" ‘’j 


i + 1 


A',  .sample  .\, 

J i 


for 


f rom 


from  BNP(N,  p)  and 


r 


13. 

If 

m V 

5 go  to  15. 

14. 

If 

^2 

satisfies  A"  , 

sample  X^ 

go 

to  21 

• 

15. 

A ♦ 

19. 

16. 

r 

17. 

Jll 

Ji+1  » 1, 

. . . , k-1. 

18. 

Jk 

^ r. 

19. 

P 

/ (1  - q). 

20. 

Sample 

X.  from  BI(N, 

P). 

21. 

N -• 

N - 

■‘j.  ■ 

22. 

If 

N = 

0,  deliver  X. 

23. 

If 

i = 

k - 1,  X.  -e  N 

J k 

and  deliver 

24. 

q 

q + 

'’Ji- 

25. 

i - 

i + 

1. 

26. 

P 

'’Ji 

/ (1  - q). 

27. 

Go 

to  A. 

The  reader  will  note  the  less  than  complete  specification  of 
tests  (13)  and  (21^)  in  steps  10  and  14  respectively.  This 
deliberate.  One  can  use  these  tests  or  slightly  more  conservative 
that  are  considerably  more  computationally  efficient. 

For  the  normality  criterion  in  (131,  let 

m “ i^pq  Cj^  = (q  - 

(31) 

Cj  • 0.053  + 0.027(q  - p)  + 0.054(q  - p)^  C3  = 1.3  . 


is 

tests 


v'dusi^Ut  1 hf  I'Oiui  i t iDii 


( t.' ) 

S : ' ■ 

* 

1 t"’ 

m 

wiit*  L 

•'■1 

1 mp  1 y 

(H 

) rowrittfii  111  [t*nns  ot  (U  ).  Ti\in 

t f niis 

it- ads  tt' 

' 

3 

c.-.’  - 

11  - 

11  . 

ra 

m 

S iiif  f 

lu-'  -•  dS 

in 

tU'df  f 

ti>r  (M)  ti.'  iipplv,  .uui  q - p is  L, 

; li.u 

.J  • 0 it 

2-> 


A:. sum  ini’,  that  tlio  usor  siH’i'itioit 


a'  sat  isl  it's  I'nt'  oau 


\ \\  ; hf  I’liu  1 va  ! f r.t  ifrm 


Sin-f  n 


S,  t'nf  li.is 


i 3ti 


,ni  • 1.-7. 5 . o.OOOSSi 


i.’liii  ti  If. ids  tv’  t hf  sli>;htly  mtTf  f ansf  rv.it  ivt'  cfiulitifn 


C. ; > S’: 


> . * ' 


(■ 


Npq 


- O.OOOS^'J 


’Y 


J Npti 


In  shfiL,  usliip,  (.'7^  i»  pl.uf  ol  (jj)  in  stt’p  10  ot  .ili’a't 


^;p^'^;  i np. 


it  foil  t'ws 


wr  i t f (.  1.' ' 


it  lull  M 1 
iNptl  anti 


and  ff  St  r if  t inp. 


> I 

. \ 


O.OOSi^li  .iVt'ivl  t hf  lU'fti  t.'  ft'mputf  m 


} 


i 

( 

I 

I 


) 


I 

i 


24 


-1.5m  . . , 

e . in  summary,  one  can  avoid  exponentiation  and  square  root  transfor- 

mations as  follows  for  A'  s 0.005913: 


a. 

If 

Npq 

< 25, 

BN 

does 

not  apply. 

b. 

If 

Npq 

> 25 

and 

(37) 

holds,  use  BN. 

A more  conservative,  but  more  computationally  efficient,  test  than 
the  criterion  (22)  Implies  Is  also  possible  for  approximating  the 
Poisson  distribution  by  the  normal.  Let 


1 + i_V-2h 

2h; 

8l  — g2  = 0.0544  g3  = 0.0108 

b/Tn 

g^  = 0.2743  g5  = 0.0065  . 

For  specified  A"  one  can  write  the  criterion  of  Interest  as 

(39)  4"  2 i,  . !i+ 55  + ,(h)  . 

' h h2  h5 

If  one  adopts  the  lower  bound  A'  ^ 0.005913  In  (34)  for  A"  as  well, 
then  h S 12.1,  or  equivalently  p > (12.1)^  = 146.71,  guarantees  (39). 


Y(h) 


■( 


Consider 


(40) 


h*  = min 


A"  > 4 + ^ + Y(h) 


- 

i 


I 


I 


For  A"  > 0.05913 


h*  3.8.  When  h i h*  one  can  write  (39)  equlva- 


lent ly  as 


t I 


) 


(41) 


A"li^ 


- v(h)li^ 




2* 

+ fij*'  + 


Siiu-e  one  can  slu'w  that  >(li)h‘'  monot  on  i ea  1 1 v iKnu-asi's  t ei'  h ■ li*. 
a slightly  more  eonservat Ive  criterion  than  (41)  i« 

/ '' 

/ + R,1‘  + Sc 

(42)  u - ( i -*  5 

v»»..2 


a"k^  - g,,i'  - i;,  - 0.U8J, 


where  Y(h)h‘^  = 0.1181  lor  h = 1.8  . 


Notlco  that  (42)  and  (4*1)  di>  away  with  t lu’  noiui  tn  I'omimti'  si|iiaii' 
roots  and  to  oxpoiu>nt  iat  o. 

In  sumiiuiry,  ono  can  avoid  exponentiation  aiul  square  roiM  tr.inslor- 


mat  Ions 

as 

follows  for 

A"  ■ 

O.OOSdl  ): 

a . 

If 

M S 146.71 

use 

UI’N  . 

b. 

If 

(3.8)^  - 14 

.44  <. 

H < 146.71 

and 

(42)  holds,  use  Bl’N 

c . 

if 

p < 14.44 

and 

(45)  holds. 

uso 

BI’N. 
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